Chaotic Phase Synchronization in Bursting-neuron Models 
Driven by a Weak Periodic Force 
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We investigate the entrainment of a neuron model exhibiting a chaotic spiking-bursting behavior 
in response to a weak periodic force. This model exhibits two types of oscillations with different 
characteristic time scales, namely, long and short time scales. Several types of phase synchronization 
are observed, such as 1 : 1 phase locking between a single spike and one period of the force and 
1 : I phase locking between the period of slow oscillation underlying bursts and I periods of the 
force. Moreover, spiking-bursting oscillations with chaotic firing patterns can be synchronized with 
the periodic force. Such a type of phase synchronization is detected from the position of a set of 
points on a unit circle, which is determined by the phase of the periodic force at each spiking time. 
We show that this detection method is effective for a system with multiple time scales. Owing to 
the existence of both the short and the long time scales, two characteristic phenomena are found 
around the transition point to chaotic phase synchronization. One phenomenon shows that the 
average time interval between successive phase slips exhibits a power-law scaling against the driving 
force strength and that the scaling exponent has an unsmooth dependence on the changes in the 
driving force strength. The other phenomenon shows that Kuramoto's order parameter before the 
transition exhibits stepwise behavior as a function of the driving force strength, contrary to the 
smooth transition in a model with a single time scale. 

PACS numbers: 05.45.Xt, 05.45.Pq,87. 19.1m 
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I. INTRODUCTION 

Since the discovery of synchronization in pendulum 
clocks by Huygens, synchronous behavior has been 
widely observed not only in physical systems but also 
in biological ones such as pacemaker cells in the heart, 
chirps of crickets, and fetal-maternal heart rate synchro- 
nization [l| . Such synchronization phenomena have been 
studied theoretically in terms of nonlinear dynamics, par- 
ticularly by exploiting oscillator models [2|. For exam- 
ple, synchronization observed in fireflies can be modeled 
using nonlinear periodic oscillators and is described as 
phase synchronization. Further, it has been indicated 
that the notion of phase synchronization can be extended 
to chaotic oscillators. This phenomenon is called chaotic 
phase synchronization (CPS) 0-(5|. 

Furthermore, synchronization phenomena in neural 
systems have also attracted considerable attention in re- 
cent years. At the macroscopic level of the brain ac- 
tivity, synchronous behavior has been observed in elec- 
troencephalograms, local field potentials, etc. These ob- 
servations raise a possibility that such neural synchro- 
nization plays an important role in brain functions such 
as perception Q as well as even in dysfunctions such as 
Parkinson's disease and epilepsy |7t4lG|. In addition, at 
the level of a single neuron, it has been observed that 
specific spiking-bursting neurons in the cat visual cortex 



contribute to the synchronous activity evoked by visual 
stimulation [Tlj : further, in animal models of Parkinson's 
disease, several types of bursting neurons are synchro- 
nized [8|. Moreover, two coupled neurons extracted from 
the central pattern generator of the stomatogastric gan- 
glion in a lobster exhibit synchronization with irregular 
spiking-bursting behavior jl2|, [l3[ . 

Hence, it is important to use mathematical models 
of neurons to examine the mechanism of neuronal syn- 
chronization with spiking-bursting behavior. As mathe- 
matical models that include such neural oscillations, the 
Chay model [l4| and the Hindmarsh-Rose (HR) model 
fl5| have been widely used. These models can generate 
both regular and chaotic bursting on the basis of slow- 
fast dynamics. The slow and fast dynamics correspond to 
slow oscillations surmounted by spikes and spikes within 
each burst, respectively. The former is related to a long 
time scale, and the latter, to a short one. 

Phase synchronization in such neuronal models is dif- 
ferent from that in ordinary chaotic systems such as the 
Rdsslcr system, owing to the fact that neuronal models 
typically exhibit multiple time scales. However, it is pos- 
sible to quantitatively analyze the neuronal models by 
simplification, for example, by reducing the number of 
phase variables to 1 by a projection of an attractor (a 
projection onto a delayed coordinate and/or a velocity 
space QilQj]])- Recently, a method called localized sets 
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technique has been proposed for detecting phase synchro- 
nization in neural networks, without explicitly defining 
the phase [US]. 

In this paper, we focus on synchronization in peri- 
odically driven single bursting neuron models, which is 
simpler than that in a network of neurons. In previ- 
ous studies, phase synchronization of such a neuron with 
a driving force has been considered both theoretically 
[2lj and experimentally In these studies, the period 
of the driving force was made close to that of the slow 
oscillation of a driven neuron. On the other hand, in 
this work, we adopt the Chay model [3| to investigate 
whether phase synchronization also occurs with the ap- 
plication of a force whose period is as short as that of the 
spikes. In particular, we focus on the effect of the slow 
mode (slow oscillation) on the synchronization of the fast 
mode (spikes). 

It should be noted that this fast driven system may be 
significant from the viewpoint of neuroscience. In fact, 
fast oscillations with local field potentials have been ob- 
served in the hippocampus and are correlated with syn- 
chronous activity at the level of a single neuron [22j . 

From intensive numerical simulations of our model, we 
find that the localized sets technique can be used to de- 
tect CPS between the spikes and the periodic driving 
force, even in the case of multiple time scales. Further- 
more, we find two characteristic properties around the 
transition point to CPS. First, the average time interval 
between successive phase slips exhibits a power-law scal- 
ing against the driving force strength. The scaling expo- 
nent undergoes an unsmooth change as the driving force 
strength is varied. Second, an order parameter R, which 
measures the degree of phase synchronization, shows a 
stepwise dependence on the driving force strength K be- 
fore the transition. That is, R(K) does not increase 
monotonically with K but includes a plateau over a range 
of K (a step), where R is almost constant. Both of these 
characteristics are attributed to the effects of the slow 
mode on the fast mode and have not been observed in a 
system with a single time scale. 

This paper is organized as follows. Section IPTl explains 
the model and describes an analysis method for spiking- 
bursting oscillations. Section Mil presents the results 
of this study. Finally, Section IIVI summarizes our re- 
sults and discusses their neuroscientific significance with 
a view to future work. 



II. BURSTING NEURON MODEL 

As an illustrative example of a bursting neuron model, 
we consider the model proposed by Chay, which is 
a Hodgkin-Huxley-typc conductance-based model ex- 



pressed as follows fl4j : 

V = gtmlh^m - V) + gkyqHVK - V) 

+9*k,cy^(Vk -V)+ g* L (V L - V), (1) 

<i = (Qoo - Q)/r q , (2) 
C = p[ml h 00 {V c -V)-k G C\. (3) 

Equation represents the dynamics of the membrane 
potential V, where Vi, Vk, and Vl are the reversal po- 
tentials for mixed Na + and Ca 2+ ions, K + ions, and the 
leakage current, respectively. The concentration of the 
intracellular Ca 2+ ions divided by its dissociation con- 
stant from the receptor is denoted by C. The maximal 
conductances divided by the membrane capacitance are 
denoted by g*, g^ v , j^ c , and <j£, where subscripts (I), 
(K, V), (K, C), and (L) refer to the voltage-sensitive 
mixed ion channel, voltage-sensitive K + channel, Ca 2+ - 
sensitivc K + channel, and leakage current, respectively. 
Finally, and are the probabilities of activation 
and inactivation of the mixed channel, respectively. 

In Eq. ([2]), the dynamical variable q denotes the proba- 
bility of opening the voltage-sensitive K + -channcl, where 
T q is the relaxation time (in seconds), and is the 
steady-state value of q. 

It should be noted that, on the basis of the formulation 
in 01 1 the variables TOoo, ^ooj and goo are described by 
Uoo = oty/ (a y + /3 V ), where y stands for to, h, or q with 

a m = 0.1(25 + V)/[l-exp(-0.1V- 2.5)], 

p m = 4exp[-(F + 50)/18], 

a h = 0.07exp(-0.(W-2.5), 

p h = l/[l + exp(-0.iy-2)], 

a q = 0.01(20 + V)/[l-exp(-0.1F- 2)], 

P q = 0.125 exp[-(V + 30)/80]. 

Further, r q is defined as 

T q = [230(a g + Aj)] -1 . 

In Eq. kc, p, and Vq are the efflux rate constant 
of the intracellular Ca 2+ ions, a proportionality constant, 
and the reversal potential for Ca 2+ ions, respectively. 

In this study, a sinusoidal driving force with amplitude 
K and frequency / is added in Eq. ([1]) as follows: 

V = gfn&MVS -v) + gk,W{v K - V) 
+ gkc T ^-c(VK-v) + gi(v L ~v) 

+ K sin2n ft. (V) 

In the following sections, we fix the frequency at / = fo 
and vary the amplitude K to investigate the response of 
the system. 

The values of the reversal potentials and the fixed pa- 
rameters used in our simulation are listed in Table HI The 
value of g^ c significantly influences the dynamics of the 
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system, and a chaotically bursting behavior can be ob- 
served in the vicinity of g^ c = 11. We use this value 
in our simulation. Figure Q] shows (a) the chaotic attrac- 
tor, (b) the time series of V(t), and (c) the average power 
spectrum of the time series for K = 0, where the time se- 
ries includes two time scales — one for the spikes within 
each burst and the other one for the bursts themselves. 

For discussions in the rest of this paper, we introduce 
the following terminology for the time series of the Chay 
model: the fast mode describes the spiking oscillation 
in the dashed rectangles in Fig. (Hb), whereas the slow 
mode describes the oscillation in the lower envelope of 
V(t) between the dotted lines, where the dashed-dotted 
curve shows the slow oscillation. The variable C domi- 
nates the slow dynamics with the time constant 1/p. In 
fact, the decrease in V (i.e., hyperpolarization) between 
the bursts shown in Fig. [lib) corresponds to the decrease 
in C shown in Fig. HJa). On the other hand, the increase 
in C shown in Fig. Ufa) corresponds to the spiking of 
V in Fig. [TJb). Hereafter, the period for hyperpolariza- 
tion of V is called the quiescent period, as indicated in 
Fig. [Ub). It should be noted that the amplitude of K 
takes small values in our simulation as compared with 
the change in voltage for a spike, such that the driving 
force is weak. 

A clear peak can be observed in the high-frequency 
part of the power spectrum (Fig. [He)), corresponding to 
the fast spiking activity. Additionally, a broadband peak 
for the slow oscillation is observed in the low-frequency 
part. In what follows, let us investigate the system under 
force with a frequency close to that of the spiking mode 
(i.e., / = /* « 0.924 Hz), which is the natural frequency 
of the fast dynamics of the system. 



TABLE I: Parameters used in the numerical simulations. 



Parameter 


Value 


Unit 


Vk 


-75 


mV 


Vi 


100 


mV 


Vl 


-40 


mV 


Vc 


100 


mV 


V n 


-30 


mV 


Vm 


-50 


mV 


g*K,v 


1700 


s- 1 


9*i 


1800 


s- 1 


9l 


7 


s' 1 


kc 


3.3/18 


mV 


P 


0.27 





Next, in order to investigate phase synchronization be- 
tween the spikes and the periodic driving force, we con- 
sider the phase n of the driving force at t = t n when 
the nth spike occurs in the neuron (see [23j]). Figure [5] 
illustrates the manner in which the phase variable of the 
driving force at t n , which is defined as the moment when 
V(t) exceeds a certain threshold Vth, can be measured. 
Once a sequence of t n is determined, we can assign points 
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FIG. 1: (a) Chaotic attractor in the (V^, C) space. The arrows 
indicate the direction of the trajectory, (b) Time series of 
V(t) for the chaotically bursting behavior described by Eqs. 
(l)-(3). The terms "fast mode" and "slow mode" describe the 
spiking oscillation and the oscillation in the lower envelope of 
V(t) between the dotted lines (shown by the dashed-dotted 
curve), respectively, (c) Power spectrum for time series of 
V(t). The spectrum is averaged over 100 time series with a 
length of 2 16 x At, where At = 0.005 s is the time step for 
the numerical integration. 



8 n on a unit circle, where each n is determined as the 
phase of the sinusoidal force at time t n . We assume that 
< 9 n < 27r. We term these points the spiking time 
points (STPs). 

We can detect the CPS between the driven system and 
the driving force as a localization of the STP distribution; 
that is, there is an open interval on the unit circle where 
no spiking time point is detected. In [l8|, such a 
localization of the STPs (obtained from a sufficiently long 
time series) is mathematically described in the following 
manner. Let the distribution of the STPs be included in 
a set T> on the unit circle; T> is localized if there exist 
open sets A^ on the circle such that V n J2i ^ = ^ 
should be noted that only one parameter Vth is required 
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to detect the CPS by this algorithm. 
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FIG. 2: (Color online) Schematic illustration of spiking time 
points (STPs). Under the condition of upward crossing, the 
relevant threshold Vth for the spikes (bottom left) is indicated 
by the horizontal dotted line. In crossing the threshold, the 
numbered crosses (x) correspond to the plus symbols (+) on 
the sinusoidal driving force (top left). The plus symbols on 
the circle (right) correspond to those on the top left. The plus 
symbols on the circle are the phases of the STPs. 
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FIG. 3: (Color online) Time series of V(t) (red solid line) for 
driving forces (blue dotted line) with different amplitudes: (a) 
K = 0.01, (b) K = 0.113, (c) K = 0.2. The corresponding 
STPs appear on the unit circle in (d), (e), and (f), respec- 
tively. The length of the time series for plotting the STPs is 
1000 s. 



III. RESULTS 



In this section, we mainly show the results for the pe- 
riodic driving force with the frequency at /o = 0.9 Hz. 
This value is close to the natural frequency /* of the 
fast mode. In this case, the quiescent periods of V dis- 
appear when the forcing amplitude K increases; that is, 
V exhibits spiking without quiescent periods. Then, we 
find the CPS between the single spikes and the driving 
force. Moreover, we observe two characteristic phenom- 
ena around the transition point to the CPS. One phe- 
nomenon shows a power- law scaling against K, which 
is exhibited in the average time interval between suc- 
cessive phase slips. The scaling exponent undergoes an 
unsmooth change as K is varied. The other phenomenon 
shows a stepwise behavior observed in Kuramoto's or- 
der parameter for the driving force strength before the 
transition (shown in Fig. [5]). 

In addition to the results of the case / = /o, we also 
show brief results for the periodic driving force with fre- 
quency / = /i = 0.92 Hz. In this case, the quiescent 
periods of V do not disappear (shown in Fig. II 1 [) , even if 
the value of K increases to its maximum value of K « 0.5 
for the case / = /o- As will be explained later in de- 
tail, a stepwise behavior can be observed as well by us- 
ing another observation variable. All the characteristic 
phenomena are inherent to systems with multiple time 
scales. 



A. Detection of Phase Synchronization 

Detection by Spiking Time Points 

Figurcs[3|a)-[3tc) show the time series of V(t), together 
with the sinusoidal driving force. For a weak driving 
force with K = 0.01, as shown in Fig. G3a), there is no 
synchronization between the spikes and the force. On the 
other hand, when K = 0.113, as shown in Fig. [3fb) , the 
system exhibits phase synchronization, i.e., a one-to-one 
correspondence between a single spike and one period of 
the driving force. Additionally, there are fluctuations in 
the inter-spike intervals (ISIs). Therefore, this state is 
considered to be one that exhibits CPS. In the case of 
a strong driving force with K = 0.2, as shown in Fig. 
[3fc), classic phase locking (CPL) is observed with two 
periodically alternating ISIs. More precisely, each pair 
of successive spikes is observed at the same points in the 
period of the force. 

While the state shown in Fig. [3f a) does not exhibit 
phase synchronization in terms of STPs, the states shown 
in Figs. |3|b) and 3(c) exhibit phase synchronization. 
Figures [3Jd)-3(f) show the STPs on the unit circle for a 
certain time interval. The length of the time series for 
plotting the STPs is 1000 s. This length of the time series 
is sufficiently long to determine the localization of the 
distribution of the STPs, as mentioned in the definition 
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of the localization of tho STPs. In Fig. [3](e) , when K = 
0.113, we can detect CPS between the spikes and the 
force, because the STPs are localized yet distributed on 
the unit circle. As shown in Fig. [3Jf), when K = 0.2, the 
CPL can be detected on the basis of the fact that all the 
STPs are concentrated at two points on the unit circle. 
This means that each pair of successive spikes completely 
synchronizes with two periods of the force. However, 
no synchronization can be detected in Fig. [3jti) , when 
K = 0.01, because the STPs are not localized. In other 
words the entire circle is filled with STPs. 
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FIG. 4: (Color online) (a) Time evolution of the phase differ- 
ence A9 for K = 0.01 (blue dashed line), 0.064 (green dotted 
line), 0.073 (red solid line), (b) Time series of V(t) (red solid 
line) with the periodic force (blue dotted line) for K = 0.064. 
Phase slips are observed in the region indicated by the dashed 



Detection by Phase Difference 

To confirm the CPS, as shown in Fig. Elb), from an- 
other perspective, let us define the phase difference AO 
between the spiking oscillation and the driving force and 
then observe its time evolution for different values of K 
around the transition point. Suppose that at each in- 
stance when the membrane potential exceeds the thresh- 
old value, the phase variable of the spiking oscillation, 
8 S , increases by 2ir. The instantaneous phase variable of 
the external force, e , is determined at the spiking time 
without taking it to be modulo 27r unlike 9 n considered 



in the case of STPs in Section [TTJ Thus, the phase dif- 
ference is defined as AO = C — S . It should be noted 
that I : m phase synchronization between spikes and 
the external force can be defined as \mO s — 10 e \ < const. 

Figure HJa) shows the time evolution of AO for K = 
0.01, 0.064, and 0.073. For K = 0.01, the time evolution 
of AO shows an oscillation with decreasing tendency. For 
K = 0.064, A9 temporarily fluctuates within a bounded 
region (plateau) but sometimes exhibits a sudden phase 
slip. Finally, for K = 0.073, Ad always fluctuates within 
a bounded region; that is, phase slip does not occur. It 
should be noted that there exists a transition point to 
CPS near K = 0.073. Therefore, we can confirm that the 
state shown in Fig. [3Jb) represents CPS in the sense of 
the conventional definition, as well, given that K = 0.113 
is beyond the transition point. It should also be noted 
that phase slips occur when the quiescent period of V 
takes a relatively long time. These slips are indicated by 
the dashed arrows in Fig. |U[b). 
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FIG. 5: Bifurcation diagram of 8 n with respect to K. CPS 
and CPL represent chaotic phase synchronization and classic 
phase locking, respectively. 



B. Route to Phase Synchronization and Phase 
Locking 

Let us clarify the route to phase synchronization by 
the dynamics of n for the driving force at the nth spike, 
which is an angle of the STPs on the unit circle. The 
bifurcation diagram for n with respect to K is shown 
in Fig. O For values of K that are relatively small, n 
can take any value in the range between and 2ir, which 
means that CPS does not occur (a non-CPS state). After 
the first transition at K = K c « 0.073, the value of 9 n is 
confined within a localized region. 

First, wc examine the system behavior around K = 
K c . Figure [S] shows the return plots in the space 
(6 n , 6 n+ i). The points shown in Fig. H^b) (K = 0.0735) 
are distributed within a limited region, whereas the 
points shown in Fig. Uta) (K = 0.069) are distributed 
throughout the entire space. Therefore, these plots im- 
ply that the system undergoes a boundary crisis 24 1. 
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FIG. 6: Return plots of 8 n in the space (8 n , 8 n +i) for (a) 
K = 0.069, (b) K = 0.0735, (c) K = 0.11, and (d) K = 0.115. 



Here, the crowding of points shown in Fig. EJb) repre- 
sents CPS (a CPS state) as in the case of the localization 
of the STPs on the unit circle, implying that phase slip 
does not occur. 

Moreover, the second transition at K s=s 0.11205 seems 
to be an interior crisis pij . which is indicated by the 
change between Fig. [6jc) (K = 0.11) after the crisis and 
Fig. EJd) (K = 0.115) before the crisis. More precisely, 
the two disconnected attracting sets shown in Fig. Iljd) 
are included in the single attractor shown in Fig. Etc). 
After the second transition, a typical sequence of inverse 
period-doubling bifurcations occurs. As a result, CPL (a 
CPL state) is observed in the regime where the driven 
system fires periodically. 
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FIG. 7: (a) Average time intervals between phase slips with 
respect to K around the phase transition point, (b) Log- log 
plot of the average time intervals with respect to the difference 
between the parameter K and its critical value K c ~ 0.073. 
The slopes represent the scaling exponents. The results are 
averaged over 100 and 1000 different realizations for (a) and 
(b), respectively. 



C. Characteristics around the Transition Point 



Scaling Behavior for Phase Slips 

We characterize the average time interval (r s ; ) between 
two successive phase slips (i.e., the plateau length) with 
respect to K around the transition point, as shown in Fig. 
[Zja) J25|. Figure[7£b) shows the log- log plot of (t s i) in de- 
pendence on \K—K C \, where K c « 0.073 is the transition 
point to CPS. Here, K c is approximately determined by 
the point at which the distribution of the STPs begins to 
become localized. We numerically find a power-law scal- 
ing (t s i) ~ \K — A' c | 7 , where the scaling constant sud- 
denly changes from 7a — 1.5 to 7 « — 7.5. Although the 
range of the scaling region is relatively short, this scaling 
behavior is clearly different from that observed for a sys- 
tem with a single time scale. The scaling behavior with 
a single time scale is related to type-I intermittency and 
is described by log(r si ) - \K ~ A c | -1 / 2 [HH^. 

In general, for a periodically driven chaotic system 



with a single time scale and a single rotation center, a 
simplified mapping model can explain the transition to 
CPS between the system and the driving force. In the 
map model, the boundary between the synchronization 
state and the non-synchronization state is explained by 
a saddle-node bifurcation of unstable and stable periodic 
orbits of the map 23|. However, in the present system, 
phase locking cannot simply be related to a saddle-node 
bifurcation, because multiple time scales exist. That is, 
the spiking period is characterized by fast oscillations, re- 
lated to the variables V and q, and the quiescent period is 
characterized by slow oscillations, related to the variable 
C. Hence, the dynamics of the system on the thresh- 
old, which corresponds to the map model, is affected by 
both the fast and the slow oscillations. Therefore, the 
mechanism of a sudden change in the scaling law differs 
from the case of a single time scale but is still an open 
problem. 
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FIG. 8: (a) Order parameter R as a function of K. The 
dashed line indicates K c or the transition point to CPS. (b) 
Enlargement of (a) for < K < 0.1. It should be noted that 
F denotes the step region before the transition to CPS. The 
results are averaged over 1000 different realizations. 



Stepwise Behavior Measured by Order Parameter 

The distribution of the STPs is now used to detect 
phase synchronization. In order to quantify the degree 
of phase synchronization, we employ Kuramoto's order 
parameter defined as R — \ Y^j=i cxp(i6 , J - ) /A^| , where N 
is the number of STPs [27j |. The parameter R satisfies 

< R < 1. This corresponds to the amplitude of an 
average of STPs. The more localized the distribution of 
the STPs is, the greater is the value of R. 

Figure EJa) shows the order parameter R as a function 
of K, with R averaged over 1000 different realizations for 
a time interval of 1000 s. The transition point K c to CPS 
is denoted by the dashed line in Fig. [HJa). We find that 
a stepwise behavior precedes the transition to CPS, with 
the step region T, as shown in Fig. |5Jb). This stepwise 
behavior indicates that there exists a region between the 
non-CPS state and the CPS state, where the degree of 
synchronization is not sensitive to K. In addition, R 
decreases just before the transition point. To the best of 
our knowledge, such a stepwise behavior in the transition 
has not been observed for CPS in coupled systems with 
a single time scale. In what follows, we will explain how 
the stepwise transition and the decrease in R are related 
to the existence of both the slow and the fast dynamics in 
the present system. In particular, we will investigate the 
probability density distribution of 6 n at K = K^, with 

1 = 1, . . . , 5, as indicated in Fig. [5th). 
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FIG. 9: (Color online) Probability density distributions of 9 n 
for K = 0.01 (red solid line) and K — 0.05 (green dashed 
line). 



Figure |9] shows the shape of the probability density 
distributions P(9 n ) of the STPs on the unit circle, i.e., 
the distributions of 9 n for K = K2 = 0.05 and for K = 
K\ = 0.01. It should be noted that CPS is not detected 
in either case; that is, the points are well distributed 
around the circle. However, we can find a peak in the 
probability density distribution for K = K% , whereas the 
distribution for K = K\ is almost uniform. 

The appearance of the peak, as shown in Fig. [9j and 
a shift in the peak in the step region T, as explained in 
Appendix \K\ are two consecutive stages that constitute 
the entire stepwise phenomenon. In the first stage, with 
an increase in K from zero to a value near K2, a peak 
in P(0 n ) of the STPs emerges near 9 n = 0.5. Then, in 
the second stage, the position of the peak shifts because 
of the effect of slow dynamics, i.e., an increase in the 
number of short inter-burst intervals. During the shift of 
the peak, R does not change very significantly because 
the value of K varies in T. Thus, the step region can be 
observed. These two stages are investigated in detail in 
Appendix [A] 



D. CPS with Quiescent Period 

In the above results, we have primarily investigated the 
1 : 1 phase synchronization between the spikes and each 
period of the driving force, wherein quiescent periods dis- 
appear for the values of K after phase synchronization. 
On the other hand, when the frequency of the driving 
force is / = /1 = 0.92 Hz, the quiescent periods do not 
disappear even if the amplitude K increases to the same 
level as in the case / = /o = 0.9 Hz. It should be noted 
that the 1 : 1 CPL can be observed for a sufficiently 
large K, e.g., K = 2. For values of K < 0.5, 1 : 1 phase 
synchronization cannot be clearly observed between the 
spikes and the driving force. However, if we define 9„ 
at the time when V approaches the minimum voltage in 
each specific quiescent period, where V decreases under 
the threshold V* = —47, we can detect CPS in the sense 
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FIG. 10: (a) Bifurcation diagram of n with respect to K. 
(b) The corresponding values of the order parameter R. The 
inset shows a magnified image for < K < 0.1. The results 
are averaged over 1000 different realizations. 



of localization of the distribution of 0„ on the unit cir- 
cle. It is important to note that the change in the order 
parameter (derived from 0„) with respect to K does not 
depend on the value of V*, sensitively. In other words, 
the value of V* that is less than approximately —46.4 
causes a stepwise transition in the order parameter for 
/o = 0.92 Hz using 0„, as shown below. 

Figure ITUTa) shows the bifurcation diagram of 8 n de- 
pending on K. We also observe a stepwise behavior be- 
fore a transition to the CPS in the variation of the order 
parameter computed from 0„, as shown in Fig. HOf b). 
Moreover, in the region of CPS, fine tuning of K yields 
1 : I CPL between a burst and I periods of the force, 
where I = 14, 15, 18, 19, etc. Figure Qj] shows 1 : 19 
CPL for K = 0.302. 



IV. SUMMARY AND DISCUSSION 

We have investigated CPS in a spiking-bursting neuron 
model under periodic forcing with a small amplitude but 
with a frequency as high as that of the spikes. First, we 
observed 1 : 1 CPS between the spikes and the periodic 
force. This CPS has been detected on the basis of the 
fact that a set of points, which are conditioned by the 
phases of the periodic force at each spiking time, was 
concentrated on a sector of the unit circle. 
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FIG. 11: (Color online) 1 : 19 CPL between bursting (red 
solid line) and the periodic force (blue dotted line). 



In addition to CPS, we observed two characteristic 
phenomena around the transition point to CPS. One phe- 
nomenon involves a change in the power-law scaling for 
the average time intervals between phase slips, as shown 
in Fig. EJb). This scaling behavior is different from that 
exhibited by the conventional system, i.e., a chaotic sys- 
tem whose attractor has a single rotation center with only 
one characteristic time scale. In such a conventional sys- 
tem, the scaling exponent for the transition to CPS takes 
a unique value of —1/2. This might be because the phase 
synchronization in question cannot be simply associated 
with a saddle-node bifurcation, owing to the interaction 
between the slow and the fast dynamics. The other phe- 
nomenon shows a stepwise behavior before the transition 
to CPS, as shown in Fig. [Hfb) . This phenomenon has 
been found by the observation that the degree of phase 
synchronization is not sensitive to the amplitude of the 
force just before the transition point. Moreover, we found 
that a decrease in the degree of synchronization appears 
(at K = if 4 in Fig. [5Jb)) even if the amplitude of the 
force is increased. The stepwise behavior and this de- 
crease could be induced by the effect of slow dynamics 
(see Appendix lA} . 

From the viewpoint of ncuroscicncc, our system might 
be regarded as a simple model whose fast driving force 
corresponds to sharp wave-ripples. This phenomenon 
involves a very fast oscillation of local field potentials 
observed in the hippocampus in the brain [22[. Let us 
discuss, below, a possible interpretation of our results 
in terms of synchronization phenomena in real neuronal 
systems. 

First, the observed phenomena in our model, particu- 
larly CPS with a fast driving force, can be interpreted as 
a consequence of the interaction between the ripples and 
the activity of a single neuron. In fact, it has been exper- 
imentally shown that such ripples synchronize in phase 
with the spikes of a single neuron [28l - t30l | . 

Furthermore, it has been suggested that firing se- 
quences accompanying ripples in the hippocampal net- 
work form a representation of stored information. The 
replay of the firing sequences during sleep mediates a 
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consolidation of memory for the stored information in 
the hippocampus and the neocortex [3iT[34l ] . In addi- 
tion, some experimental results have indicated that such 
a memory replay is conducted on a shorter time scale 
than the actual experience, where the spatiotemporal 
structure of the firing sequences plays a key role in the 
stored information [351438 1. 

Thus, if our spiking-bursting system with a fast driv- 
ing force can be regarded as a model for such biological 
neurons, our result would imply that the precision of the 
temporal structure of the spiking patterns might be en- 
hanced by CPS in real neuronal systems. It should be 
noticed that CPS is not affected even if there exist small 
fluctuations in ISIs and that CPS extends the detection 
of the temporal structures of the firing patterns in the 
short time scale of a spiking-bursting behavior. 

Additionally, from another point of view, the slow os- 
cillations along the lower envelope of the membrane po- 
tential V can be considered as transitions between the 
UP and DOWN states in a cortical neuron. Specifically, 
we focus on the UP and DOWN states during the slow- 
wave sleep. Here, the neurons in the UP state fire syn- 
chronously with higher frequency, whereas the activity 
of the neurons in the DOWN state is relatively quies- 
cent [39j . From this perspective, the disappearance of 
the quiescent periods (DOWN states) in the forced spik- 
ing in our simulation may be interpreted as a persistently 
depolarized UP state observed for cortical neurons in the 
awake states. In fact, the prolonged DOWN states during 
sleep are likely to occur owing to a decrease in the excita- 
tory input |40j. Therefore, if this input can be regarded 
as our sinusoidal driving force, CPS may efficiently help a 
weak input to depolarize a DOWN state into a persistent 
UP one. Recently, Ngo et al. reported that a neural net- 
work model can qualitatively reproduce the experimental 
results in [39j using a time-discrete map model, which is 
simpler than our model j4l| . 

Therefore, we may infer that our results on CPS at high 
frequency in the simple spiking-busting neuron model 
provide some suggestions for neuroscience to understand 
the mechanisms of the abovementioned real neuronal ac- 
tivity, particularly in terms of nonlinear dynamics. 

The following topics might be considered from the 
viewpoint of extending this study. First, for the future 
study it would be important to confirm whether our find- 
ings can be observed in other slow-fast models such as the 
Hindmarsh-Rose model fl5| and other biophysical mod- 
els (42M44| . A comparison with the other models will 
provide further insight into slow-fast dynamics in the 
neuronal spiking-bursting activity. Moreover, it should 
be important to clarify the general onset mechanism of 
the observed phenomena using map-based models, as was 
clarified in [23|, IH| • It should also be important to inves- 
tigate the response of coupled bursting systems to sinu- 
soidal forcing in terms of the interaction between the slow 
and the fast dynamics. 



Acknowledgments 

The authors are grateful to Prof. M. Tatsuno, Prof. H. 
Hata, Prof. M. Baptista, Dr. K. Morita, Dr. S. Kang, 
and Dr. G. Tanaka for their valuable suggestions. This 
work is partly supported by Aihara Complexity Mod- 
elling Project, ERATO, JST, and the Ministry of Educa- 
tion, Science, Sports and Culture, Grant-in- Aid for Scien- 
tific Research No. 21800089 and No. 20246026, and the 
Aihara Project, the FIRST program from JSPS, initiated 
by CSTP, and grants of the German Research Foundation 
(DFG) in the Research Group FOR 868 Computational 
Modeling of Behavioral, Cognitive, and Neural Dynam- 
ics. 



Appendix A: Mechanism of Stepwise Behavior 

In what follows, we explain how R increases toward 
the step region in Fig. [5Jb) and how it does not increase 
in the step region T. These phenomena can be explained 
on the basis of the changes in peaks in the probability 
distribution P(9 n ) as follows. 

First, we investigate the quiescent periods in terms of 
the extent in the decrease in V, i.e., the extent of hy- 
pcrpolarization, with an increase in K. As shown in 
Fig- H2f a). we observe two types of hypcrpolarizations 
in the quiescent periods, namely, shallow and deep hy- 
pcrpolarizations, when Uw — 46. Shallow hypcrpolariza- 
tions are observed between V = —46 and —46.5, whereas 
deep hyperpolarizations are observed below V = —46.5. 
Therefore, the two types of hyperpolarizations are dis- 
tinguished by two appropriate thresholds for V, i.e., 
Vthi = —46 and Vtya = —46.5, as shown in Fig. IT27 a). 
The threshold Vthi detects both the shallow and the deep 
hypcrpolarizations, whereas the threshold Vthi detects 
only deep hypcrpolarization. Figure fT2T b) shows the ra- 
tio of shallow and deep hyperpolarizations over all the 
hyperpolarizations detected by Vthi as a function of K, 
when counted in the interval of 1000 s and summed over 
100 different realizations. 



R on the increase toward the step region V 

As shown in Fig. I12f b). the number of deep hypcrpo- 
larizations (blue dashed line) decreases with an increase 
in K. This implies that the number of longer ISIs corre- 
sponding to deep hypcrpolarizations decreases, and this 
decrease is related to an emergence of a peak in P(6 n ), 
as explained below. 

Let x n be the length of the ISI between the nth spike 
and the [n + l)th spike. Then, Figs. IT37 a)((c)) and 
(b)((d)) show the time series of 9 n and x n for K = 
(K = 0.05), respectively. In Fig. [T37c). we observe grad- 
ually shifting envelopes for the oscillations of 9 n , as indi- 
cated by the dotted rectangles, whereas no such specific 
tendency is observed in Fig. [TBT a) . The gradual shifts in 
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FIG. 12: (Color online) (a) Time series of V(t) with two 
thresholds, Vthi = —46 (blue dashed line) and Vth2 = —46.5 
(green dashed-dotted line), for K = 0. (b) The ratio of shal- 
low (green dashed-dotted line) and deep (blue dashed line) 
hyperpolarizations with respect to K. The results for (b) 
are averaged over 100 different realizations. Also shown is R 
(red solid line). 
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FIG. 13: Time series of n (a)((c)) and x n (b)((d)) for K = 
(K = 0.05). The regions indicated by the double-headed 
arrows in (d) correspond to the bursting periods other than 
the period of deep hyperpolarizations. 



the envelope occur during the time intervals other than 
those during which deep hyperpolarization is observed, 
which arc indicated by the double-headed arrows in Fig. 
ITBT d). It should be noted that deep hyperpolarization 
corresponds to the oscillations of x n exceeding x n ~ 2.5 
and that Fig. ITBT b) shows a greater number of deep hy- 
perpolarizations than those shown in (d). 

These observations give rise to a change in the shape of 
P(9 n ) in Fig. M That is, P(6 n ) is uniformly distributed 
for K = 0, owing to the fluctuation of 9 n with a few grad- 
ually shifting envelopes, as shown in Fig. [TSTa). On the 
other hand, P(9 n ) has a peak for K = 0.05, attributed to 
the gradually shifting envelopes, as shown in Fig. I13f c). 
The region for the envelopes is around 9„ m 1 and corre- 
sponds to the peak of P(9 n ) in Fig. [§] This appearance 
of the peak in P(9 n ) is attributed to the decrease in the 
number of deep hyperpolarizations, as seen in Figs. [T3lb) 
and (d). 

R on a plateau in the step region T 

We explain how R as a function of K remains flat in V. 
Figure IT2T b) shows that the number of shallow hyperpo- 
larizations increases with an increase in K in T, as well 
(green dashed-dotted line). In order to investigate the 
increase in shallow hyperpolarizations in T, we consider 



the distribution of the ISIs. Figure Eta) shows the dis- 
tribution of ISIs for K = 0.055, 0.063, and 0.068, denoted 
by K^,K4, and K§, respectively. In this figure, we can 
observe peaks located at x„G [1.5; 2], as indicated by (i) 
in Fig. H4T a). These peaks correspond to shallow hyper- 
polarizations and become prominent at the right edge of 

r. 

Furthermore, for each K, the distribution of ISIs has 
two other peaks at smaller values of x n , as indicated by 
(ii) in Fig. [MT a). These peaks can be inferred from the 
time series of voltage V(t) for a value of K in T, as shown 
in Fig. [T31 As can be seen in the interval of t between 
350 s and 380 s, there exists a transient 1 : 3 phase 
synchronization between one burst and three periods of 
the force. The peaks located near x n G [1.5; 2] in Fig. 
I14f a) correspond to the inter-burst intervals during such 
a transient phase synchronization. In contrast, the two 
other peaks at small x n indicated by (ii) in Fig. [T4l a) 
correspond to the two ISIs for the three consecutive spikes 
in each burst. 

We also calculate the corresponding distributions of 9 n , 
as shown in Fig. HH b). The peaks in the distribution 
of related to the peaks in P(9 n ). The second 

and the third spikes in each burst during the transient 
1 : 3 phase synchronization correspond to the peaks 
in P(9 n ) around 9 n 2 and 9 n w 5 ((iii) and (iv) in 
Fig. rniTb)). respectively. It should be noted that for 
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FIG. 14: (Color online) (a) Probability density distribution 
of ISIs normalized with respect to their length x„ for K — 
0.055 (= K 3 ) (red dashed line), K = 0.063 (= Ka) (blue solid 
line), and K — 0.068 (= K$) (green dashed-dotted line), (b) 
Probability density distribution of 8„ corresponding to (a). 



the amplitude of the average of the STPs whose distribu- 
tion on the unit circle reflects the shape of P{9 n ). There- 
fore, the value of the order parameter R does not increase 
in r. However, as K increases further, the peak located 
at 9 n ps 1 becomes higher, while R increases steeply, and 
the transition to CPS occurs. 

Finally, looking at Figs. IT47 a) and ll"47 b) in more de- 
tail, we find that the peaks in both P(x n ) and P(9 n ) for 
K4 are sharper than those for K3 and K$. The sharpness 
of the peaks corresponds to the duration of the transient 
1 : 3 phase synchronization. In fact, the contribution 
to R from the newly emerged peaks in P(9 n ) around 
9 n ps 2 and 9 n ps 5 is the greatest when K ps K4. On 
the other hand, the contribution of the first concentrated 
peak around 9 ps 0.5 decreases at that value of K. As 
a consequence, even though K increases, which would 
normally cause the degree of synchronization to increase, 
the value of R decreases near K = K4. A similar phe- 
.£=0.063 




360 r , 380 
t W 



FIG. 15: (Color online) Time series of V(t) (red solid line) 
with the periodic forcing (blue dotted line) for K = 0.063. 



K = A'3, the distribution of P(9 n ) is concentrated, with 
a large peak (unimodal) located near n £ [0.5; 1]. Then, 
for K = K4, that peak is weakened and other peaks 
(multimodal) appear near 9 n ps 2 and 9 n ps 5. 

The shift in the peak in the distribution of P(9 n ) from 
unimodal to multimodal does not significantly influence 



nomenon has been reported as anomalous phase synchro- 
nization, whereby coupling among interacting oscillator 
systems increases the natural frequency disorder before 
synchronization (4o| . 
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